Experimental realization of quantum algorithm for solving linear systems of equations 
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Quantum computers have the potential of solving certain problems exponentially faster than 
classical computers. Recently, Harrow, Hassidim and Lloyd proposed a quantum algorithm for 
solving linear systems of equations: given an N x N matrix A and a vector b, find the vector x that 
satisfies Ax = b. It has been shown that using the algorithm one could obtain the solution encoded 
in a quantum state \x) using 0{\ogN) quantum operations, while classical algorithms require at 
least 0{N) steps. If one is not interested in the solution x itself but certain statistical feature of 
the solution (a;|M|a;) (A/ is some quantum mechanical operator), the quantum algorithm will be 
able to achieve exponential speedup over the best classical algorithm as A'^ grows. Here we report 
a proof-of-concept experimental demonstration of the quantum algorithm using a 4-qubit nuclear 
magnetic resonance (NMR) quantum information processor. For all the three sets of experiments 
with different choices of b, we obtain the solutions with over 96% fidelity. This experiment is a first 
implementation of the algorithm. Because solving linear systems is a common problem in nearly all 
fields of science and engineering, we will also discuss the implication of our results on the potential 
of using quantum computers for solving practical linear systems. 
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The proposition of quantum computer dates back to 1980s [T], but it was not until the late 80's and early 90's 
that quantum computers are shown to be more powerful than classical computers on various specialized problems 
[2H5]. For example, the Deutsch-Jozsa algorithm [5^, Shor's quantum algorithm for factoring integers [3], Grover's 
quantum search algorithm 4 and algorithms for Hamiltonian simulation of quantum systems [5] have been found 
to require significantly less computational steps than their classical counterparts and thus render many classically 
intractable problems realistically solvable with a quantum computer. There has been experimental demonstrations of 
important quantum algorithms such as Shor's algorithm [6], Grover search algorithm [7], optimization problems [8], 
quantum simulation of molecular systems [5Hl2). all using small-scale quantum information processors. In this work 
we implement an algorithm for solving linear systems using NMR. 

Linear systems of equations play an important role in nearly all fields of science and engineering. In quantum 
reactive scattering, the Kohn variational calculation involves the inversion of the augmented stiffness matrix [TH[TS], 
which is equivalent to solving a linear system in certain occasions. In chemistry, linear equations arise commonly in 
problems such as electrostatic calculation in density functional theory, where the discretized Poisson equation takes a 
linear form |16| . Recently the Finite Element method starts to be adopted for solving electronic structure problems 
in quantum chemistry 117 1. where Schrodinger's equation is recast in form of a linear equation. Also solving linear 
systems of equations often play a role as intermediate step in many algorithms such as quantum algorithm for data 
fitting [18]. 

Here we consider the particular type of linear system where we are given an N x N s-sparse Hermitian matrix A 
with condition number k and unit vector b and we are interested not in the solution x itself but certain feature of 
X that can be written in form of x^ Mx [M is some linear quantum mechanical operator). As the size of matrix A 
grows, the size of the data sets which define the equations increase rapidly over time. For classical algorithms such as 
the Conjugate Gradient Method [20], it takes about total runtime of 0(A^SY^log(l/£)) to get the solution x when 
A is positive definite or 0{N sK,\og{l / e)) when A is not. For all classical algorithms, it is shown that the linear lower 
bound 0{N) in runtime scaling cannot be broken even we are only interested in x^ Mx rather than x itself |13j . 

The quantum algorithm proposed by Harrow, Hassidim and Lloyd [13j for solving linear systems of equations is 
shown to improve the runtime scaling to the 0(log7V) regime. The algorithm starts with a quantum state \b) and a 
few ancilla qubits in state |0). Here the state vector of |6) in the computational basis represents the unit vector b. 
Applying the well-known phase estimation subroutine [2Ij on |0 • • • 0)|6) with U — e^'"**" as the unitary operation, we 
obtain the final state of the two-register system which is approximately l3j\Xj)\uj) up to a normalization constant 
(When the eigenvalues of A can be exactly encoded using the ancilla bits, which is the case of our experiment as we 
will see later in the paper, the final state of the phase estimation is exactly proportional to Here \uj) is 

the eigenbasis of A and \b) = j3j \uj) . When the phase estimation is completed, an ancilla bit is added to the system 
and a controlled rotation is performed on it with the register as the control register. The state of the system 
then becomes ^-y/l — C^/A||0) +C/\j\l)^ l3j\Xj)\uj) where C is a normalization constant. By inverting all the 
previous quantum operations except for the controlled rotation on the ancilla bit, we uncompute the register back 
to the state |0 • • • 0). Conditioning on measuring |1) in the ancilla qubit. The algorithm probabilistically outputs a 
state \x) = I3j\j'^\uj) in the register that is initially in the state \b). The state vector of \x) in the computational 
basis is proportional to the solution x of the linear system Ax = b. The total number of quantum operations needed 
for the algorithm is 0{s^k^ logiV). This implies an exponential speedup over the Conjugate Gradient Method, which 
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is the best classical algorithm for solving the problem. 

Let us consider a 2 x 2 system to find x = {xi X2)^ such that Ax = b where 




where b is normalized. The spectral decomposition A — Y^j=i^j\^j){'^j\ has Ai = 1, = — |1)) and 

A2 = 2, \u2) = :^(|0) + The quantum algorithm to solve this problem can be summarized as the following six 
steps HBj (see Fig. [s]). 

N 

1. Input the vector 6 as a quantum state \b) — ^ bi\i) stored in a quantum register (termed register B) and prepare 

1=1 

T-l 

the register C in the state ^ k) = id^O) + |01) + |10) + |11)). Here \i) and |t) represents the computational 

T=0 

bases of registers B and C that encode the value of i with the n — log = 1 qubit and t = log T = 2 qubits 
respectively. 

T-l 

2. Apply the conditional Hamiltonian evolution ^ ® ^-^^'^^o/t both registers B and C up to error sh- 

r=0 

Here 'c' in the notation |r) {t\'^ marks the control register and tg = 2it. 

3. Apply quantum Fourier transform to the register C. Denote the Fourier basis as |fc). In general, at this stage 
in the superposition state of both registers, the amplitudes of the basis states are concentrated on k values that 
approximately satisfy « '^r^, where Afe is the fc-th eigenvalue of the matrix A. However, A is defined as in 
Eq. ([1]), the eigenvalues Ai, A2 are exactly powers of 2. Therefore if we further let = 27r, \k) — \\k), the phase 
estimation subroutine will produce exactly ft l-^j) P'^Y^'^- 

4. Add an ancilla qubit and apply conditional y-rotation Ry on it, controlled by the register C. This rotation 
transforms the qubit to ^1 — ^ \ ^) ~^ ^ 1 1) ' where (7 is a nonzero normalization constant . This is a key step of 
the algorithm and it involves finding the reciprocal of the eigenvalue \j quantum mechanically and y-rotate the 



qubit with an angle 6j such that sin6'j = To achieve both, we first use a SWAP gate (Fig. 
the transformation 



8 ) to accomplish 



|Ai> = |01) |10) - |2Ar^) 

(2) 

in register C. Then we use the register C with |2AJ^) states as the control register to apply the controlled 
Ry{\~^) rotation on the ancilla qubit. As shown in Fig. [sj the controlled- i?^ (0) gates applies rotation with 

= (27r/2'~)Aj^ where r is a parameter. In this experiment we let r = 2 and use the approximation 
sin{9j/2) w Sj/2. With this approximation, after the controlled- _Ry rotation gates, the state of the system is 
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FIG. 1. (Color online) The four-qubit quantum circuit for the algorithm. The numeric labels 1 to 6 represent the 6 steps 
outlined in the text. Here r = 2 and to = 27r, determining the precision and probabihty of the correct answer. The matrix 

I 0\ J_ f 1 1\ p _ / cos(0/2) -sin(e/2) 

0V ^/2Vl-l/ ^ V sm(6»/2) cos(e/2) 



forms of the single-qubit gates are S ■ 
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FIG. 2. (Color online) Properties of the molecular iodotrifiuoroethylene. The chemical shifts and J-coupling constants (in Hz) 
are on and below the diagonal in the table, respectively. The chemical shifts are given with respect to reference frequencies of 
376.47 MHz (fluorines) and 100.64 MHz (carbons) at 303. OK. The transmitter offsets of carbon channel and fluorine channel 
are set at 15479.7Hz and -44701.0Hz respectively. The molecule contains four weakly coupled spin half nuclei which are 
13q^ i^_Fi, ^^F2, ^^Fs. The natural abundance of the sample with a single ^^C is about 1%. 



close to 




(3) 



and C is 0.736 according to the calculation based on Ref . |22| . 
5. Uncompute register B and C by inverting the operation in the steps 1-3. 



N 



6. Measure the ancilla bit. If it returns 1, the register B of the system is in the state ^ fijX - ^ \uj) up to a 

3=0 

normalization factor, which is equal to the solution x of the linear system Ax = b. 

The experiment was carried out on a Bruker AV-400 spectrometer {9.4T) at 303. OK. We chose iodotrifiuoroethylene 
dissolved in d-chloroform, where a ^^C nucleus and three ^^Fi nuclei consist of a four-qubit quantum system. We 
label ^^C as the first qubit, ^^Fi, ^^F2 and ^^^3 as the second, third and forth qubit. Fig. 2 shows the measured 
properties of this four-qubit quantum system 23J. We control the fluctuation of the temperature to be within O.IK 
such that the effect of the chemical shift variations is suppressed. This system is first prepared into a pseudo-pure 
state (PPS) pa—^^j^I + e |0000) (0000| with / representing the 16 x 16 unity operator and e « 10^^ the polarization, 
using the line-sclective-transition method [23l |2l] • It needs two gradient ascent pulse engineering (GRAPE) pulses 
[50] and two gradient field pulses. We apply quantum state tomography [21] to reconstruct its experimental density 
matrix of the prepared PPS and the experimental fidelity is around 98.7%. The state fidelity is calculated by F = 
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FIG. 3. (Color online) Experimental spectra and reconstructed density matrices for the final states, (a), (b), (c) Experimental 
^'^C spectra of the final states after a 7r/2 readout pulse for three different fo which are listed in Fig. [4] There are eight peaks for 
carbon. Here we only show four of them related to the solution, and their intensities represent the respective probabilities of 
projecting the final state onto the states |0011) , |0000) , |0001) , |00I0) from left to right [32]. The other four peaks are almost 
zero which are not shown here. The vertical axes have arbitrary but the same units. The numbers above the peaks are the 
relative intensity compared to the intensity of the peak of PPS. The ratio of the intensity of the peaks related to |0001) and 
|0011) approximates \xi/x2\'^ of the solution x = {xi 2:2)^. We can obtain these ratios from (a), (b), (c) which are about 
1:2, 3:1, 1:1 respectively. These agree with the theoretical values in Fig. |4] The experimentally measured, fitting and ideal 
spectra are shown as the blue, red and green curves, respectively, (d), (e), (f) Real parts of experimentally reconstructed 
density matrices of the final states in the subspace where the first and the second qubits are in the |00) state, along with the 
theoretical expectations (g), (h), (i). The rows and columns represent the standard computational basis in binary order, from 
|00) to |11). The intensities of the rest parts of the real parts and all the image parts are less than 3% which can be seen in 
the supplementary material. 



Tr (pthcoiyPcxp) /y^^ (Pthcoryj (Pcxp) ' where pcxp and pthcory represent experimentally measured density matrices 
and the theoretical expectation, respectively. Then, we perform a rotation operation Ry{6) = e^*^»^ (i.e., a rotation 
along the y axis with an angle 6 to the ^^^3) [25^27) . to obtain the initial state p-m = |0006) (0006|, where the 
normalized state \b) = cos(6'/2)|0) + sin(6'/2)|l) with the state vector 6 = [cos(6'/2), sin(6'/2)]'^. 

Now, we implement the quantum circuit of the algorithm shown in Fig. [8] on the prepared input state pi„. 

The quantum circuit is realized as a whole block using a shaped radio-frequency (r.f.) pulse that is optimized by 

the gradient ascent pulse engineering (GRAPE) algorithm (SBHSU]. The GRAPE pulse is characterized by 1500 

segments, the pulse duration of 22.5 ms, and is robust to r.f. inhomogeneities, with a theoretical fidelity 0.995. 

Finally, we measure the final state to obtain the solution of the linear equations. The desired final state is \^end) — 

a|0000) + 6|0010) +c|0001) +d|0011) with \x) = c|0) + d|l), the solution x = {xi X2f = (c df /C is encoded 
^ / 

|00£C1> 

into the state of qubit 3 in the subspace labeled by |0)i|0)2|l)4. Here we consider three different 6 by preparing 
three input states pi„ with different 9. We perform a partial state tomography to get the information of c and d (see 
Method). The experimental "'^'^C spectra are shown in the left side of Fig. Ul \xi/x2f' — \c/d^ are the ratio of the 
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FIG. 4. (Color online) Experimental results of quantum algorithm for solving the linear equations Ax — b. xthcory is the 
theoretical solution while Xcxp is the experimental solution picked up from quantum state tomography. Here = 

I * I max 

^cxp^ ^theory ^ whcre sj^xp , ^l^^^^^.y 316 thc ith elements of ^theory and a^oxp respectively. The fidelity in the table refers to 

theory max 

the final state of the whole 4-qubit system, and is measured between the experimental final state and the theoretical algorithm 
outcome. 



probabilities of projecting the final state to the states |0001) and |0011), and the relative phase between xi and X2 
can be obtained by the coherence term cd* or c*d of qubit 3. Since A and b both are real, x is real and the relative 
phase is either or tt by the sign of cd* or c*d. There may exist a global phase which can not be determined in the 
solution X, but this global phase is often not important when one estimates the expectation value of some operator 
associated with x. We list all results of the three experiments in Fig. |4j The experimental errors are about 7% which 
is measured by the formula \ Axi/xi\^^^^. Furthermore, we perform the complete quantum state tomography |31| for 
the final state in the Hilbert space spanned by these four qubits which is described in supplementary materials. The 
real parts of the reconstructed density matrices in the subspace labelled by |00)i2 are shown in the right side of Fig. 
[sj with the experimental fidelities are all above 96%. 

The errors of the experiment mainly come from the inhomogeneity of magnetic fields, the imperfection of the 
GRAPE pulses and the chemical shift variations. The imperfection of the GRAPE pulses produces about 1% error to 
the final state which can cause about 2% error to the intensities of some peaks. The experimental time is about 50 ms, 
about 10% of the coherence time (T2* in Fig. 2). Hence, the decay of the signals due to the limitation of coherence 
time is an important source of errors. The chemical shift of fluorine changes by about 3 Hz when the temperature 
varies by IK [32^, so we control the fluctuation of the temperature below O.IK to reduce the error to about 0.2%. We 
obtain the intensity of the peaks by fitting the spectra using the Lorentzian curves to reduce the error by the noise. 
The total derivation between the experimental final state and the theoretical algorithm outcome is about 3%. 

The theoretical errors of algorithm concentrate on the phase estimation in the second step and the control-rotation 
operation in the forth step. In our cases, for the reason that the eigenvalues of A which we chose satisfy = ^f^, 
the phase estimation is perfect and produces no error in theory [13j . Most theoretical errors come from the control- 
rotation operation which depends on the positive integer parameter r. This error essentially depends on the quality 
of the approximation sin a Ri a, where a is the rotation angle which is integral multiple of ^ttt acted on the ancillary 
bit. The theoretical error in this step decreases as r increases when the intensities of the counterpart peaks of the 
final states are inversely proportional to |22j . Therefore, we take a balanced choice r — 2 in the experiments, which 



causes that the theoretical error is 4% contributing to 



Axj 



and the intensities of the counterpart peaks are about 
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10% of the intensity of PPS. As presented in Fig. |4j we know that the theoretical error caused by this approximation 
is one of the major errors. 

In conclusion, we experimentally demonstrate for the first time the quantum algorithm for solving linear systems of 
equations in a 4-qubit NMR system. We acquire the solutions with errors about 7% for all the three group experiments, 
which indicates fine experimental accuracy and the validity of the algorithm. This quantum experiment can be a good 
evidence that the quantum computer can help solving common and important problems. For solving linear equation 
plays important role in many classical algorithm and classical computer, this experiment is the first step to create 
and realize more quantum algorithm. It gives us the hope to broaden the application of quantum computer. 

The application of the quantum linear equation solver could be extended to a range of applications in many fields. 
For example, the ability of using the quantum algorithm to solve Poisson equation 133) would allow quantum chemists 
to speed up electrostatic calculation in density function theory. Furthermore, since the quantum algorithm could 
be used for efficiently solving linear systems of differential equations |34| . quantum computer might prove useful for 
solving the differential equations systems that arise in technical application. 

Methods 

Hamiltonian simulation. In this experiment, since the unitary evolution exp{—iAt) is a single-qubit gate, when 
realizing the exp(— lAi) operations we are not directly implementing any Hamiltonian simulation algorithms in the 
literature (such as [35 ) which are responsible for the exponential speedup of the quantum algorithm |13j . Instead, we 
decompose the exp{—iAt) gates into elementary quantum gates using Group Leader optimization algorithm |36j . 

Finding the reciprocal of the eigenvalue Xj. The inversion technique which we use in Eq. ^ to find the 
reciprocals of Xj is applicable only to some particular type of matrices A such as the one defined in Eq. ([T]), where a 
SWAP gate combined with a proper choice of r can give us a sufficiently good approximation of the amplitude C /Xj 
in the final state Eq. ([s]) of the ancilla qubit. For the general case where the i-qubit register C holds a superposition 
of states that are approximately \Xj), one can resort to techniques such as Newton iteration for finding the reciprocals 
of Xj. The quantum algorithm implementing Newton iteration has also been proven efficient |33| . Furthermore, in 
step [4] we also used the approximation dj w sin 9j . In general this approximation does not hold but one could use 
bisection method [33 to efficiently find 9j — arcsin(C'/Aj). 

Partial Tomography. The natural abundance of the sample in which just one carbon is ^^C is about 1%. 
To distinguish those molecules against the large background, we read out all three qubits via the ^'^C channel, by 
applying SWAP gates and reading out the ^^C qubit. The solution x can be obtained from the subspace labeled by 
|0)i|0)2|l)4, i.e., the information of c and d. This can be achieved by 5 readout pulses {YEEE,YEEE * SWAP12, 
YEEE*SWAPi3, YEEE*SWAPu, XEEE*SWAPi3). Here E represents the unity operator. X and Y represent, 
respectively, a 7r/2 rotation operation along x and y axis. SWAPij denotes a SWAP operation between ith and jth 
qubits The first four readout pulses are to get the probabilities of |cp and |c?p, while the last readout pulse is to 
get the relative phase of c*d. 
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SUPPLEMENTARY MATERIALS 

APPENDIX A. PULSE SEQUENCE OF THE EXPERIMENT 

Figure 5 is the pulse sequence of the experiment. The first two GRAPE pulses and two gradient field pulses were 
used to prepare a pseudo-pure state (PPS) |0G00) (00G0|. The third GRAPE pulse is to perform the quantum circuit. 
The time of the readout pulses differs from 0.4ms to 25ms. The implementation time of circuit and the readout pulse 
is about 0.05ms in total. 




20.0 



0.4 ~ 20- 



GRAPE 3 



Readout 



FIG. 5. Pulse sequence of the experiment (not to scale). The numbers that label the pulse symbols represent the time duration 
of the pulse in milliseconds. The symbols 'C and 'F' represent the channels for Carbon and Flourine atoms respectively. The 
rectangular boxes in the figure represent gradient field pulses. 



APPENDIX B. CHANGE OF CHEMICAL SHIFT CAUSED BY TEMPERATURE FLUCTUATION 

The chemical shifts are given with respect to reference frequencies of 376.47 MHz (fluorines) and 100.64 MHz 
(carbons) at 303. OK. The chemical shifts of fluorines change by about 3 Hz when the temperature varies by IK while 
the chemical shift of carbon is almost the same. The expression for the relationship between the chemical shifts of 
fluorines and the temperature is obtained by fitting the experimental data: 

w_F, = -33122.4 - 3.0(T - 303.0) 
w_F, = -42677.7 - 1.3(T - 303.0) 
wj,3 = -56445.8 + 1.6(T - 303.0) 

where wp^, w^j, ujf-j, represent the the chemical shifts of i^i, i^2j ^3 respectively and T represents the temperature. 
The units are Hz and K. We controlled the temperature in 303.0 ± O.liiT in the experiment. 



APPENDIX C. STATE TOMOGRAPHY. 

The state tomography for the whole 4-bit quantum state needs 44 readout pulses: EEEE, EX EE, EYEE, EEXE, 
EXXE, EYXE, EEYE, EXYE, EYYE, EEEX, EXEX, EYEX, EEXX, EXXX, EYXX, EEYX, EXYX, 
EYYX, EEEY, EXEY, EYEY, EEXY, EXXY, EYXY, EEYY, EXYY, EYYY, swapi2 * EEYY, swapu * 
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FIG. 6. (Color online) Spectra and state tomography of the pseudo-pure state, (a) is the spectra of C obtained by 7r/2 readout 
pulse. The vertical axe have arbitrary unit, (b) is the real part of state tomography of PPS. (c) is the image part of tomography 
of PPS. The fidelity of the whole PPS is 98.73%. 



EEXY, swapi2 * EEEY, swapi2 * EEYX, swapi2 * EEXX, swapi2 * EEEX, swapi2 * EEYE, swapi2 * EEXE, 
swapi2 * EEEE, swapu * EEEY, swapiz * EEEX, swapi^ * EEEE, swapu * EEEE, YEEE, YEEE * swapi2, 
YEEE * swapis, YEEE * swapi^. {E represents the unity operator. X represents a rotation operation with the 
angle 7r/2 along x axis while Y along y axis, swapij means a swap operation between the ith and j'th qubits.). As 
discussed in the methods, we can pick up the accurate and complete information of x using only 5 readout pulses: 
YEEE, YEEE * swapi2, YEEE * swapis, YEEE * swapu, XEEE * swapi^ to get the solution x. Combining the 
first four readout pulses , we can obtain the diagonal elements of the density matrix. Using the fifth readout pulse, 
we can get C24 and C42. The state tomography results using 44 readout pulses for PPS and the experimental final 
states are showed in Figure 6 and Figure 7. 
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FIG. 7. (Color online) Final state tomography. (a),(c),(e) are the real parts of the state tomography of the experimental final 
states for experiment 1,2,3 respectively, while (b),(d),(f) are the real parts of the state tomography of the theoretical final 
states. The rows and columns represent the standard computational basis in binary order, from |0000) to |1111). The fidelitis 
of the experimental final states are 96. 4%, 96. 6%, 96. 7%, respectively. 

APPENDIX D. RELATIONSHIP BETWEEN PEAKS AND ELEMENTS OF DENSITY MATRIX 

As shown in Figure 8, there are eight peaks for carbon. Four of them are almost zero, and the intensity of the four 
significant peaks quantify the probabihties of the states |0001), |0000), |0011), |0010) from left to right respectively. 
In fact, defining the probabilities of the states |0000), |0001), ... , |1111) are po, pi, pi5, the four peaks which can 
be seen obviously are proportional to pi — pg, po — pg, p^ — pil, p2 ^ pio from left to right. Since pi ~ for alH > 4 
in experiment, the intensity of the four peaks are approximate proportional to the probabilities of the states |0011), 
|0000), 10001), 10010) respectively. 
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FIG. 8. (Color online)Experimental spectra of C. (a),(b),(c) are the final state spectra by a gradient field pulse and a 7r/2 
readout pulse in experiment 1,2,3 which are listed in table 2. The vertical axes have arbitrary but the same units. 



